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Abstract 

We report Monte Carlo simulations of the dynamics of a "chucker": a colloidal particle which emits 
smaller solute particles from its surface, isotropically and at a constant rate k c . We find that the 
diffusion constant of the chucker increases for small k c , as recently predicted theoretically. At large 
k c the chucker diffuses more slowly due to crowding effects. We compare our simulation results 
to those of a "point particle" Langevin dynamics scheme in which the solute concentration field 
is calculated analytically, and in which hydrodynamic effects arising from colloid-solvent surface 
interactions can be accounted for in a coarse-grained way. By simulating the dragging of a chucker, 
we obtain an estimate of its apparent mobility coefficient which violates the fluctuation-dissipation 
theorem. We also characterise the probability density profile for a chucker which sediments onto a 
surface which either repels or absorbs the solute particles, and find that the steady state distributions 
are very different in the two cases. Our simulations are inspired by the biological example of 
exopolysaccharide-producing bacteria, as well as by recent experimental, simulation and theoretical 
work on phoretic colloidal "swimmers". 
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I. INTRODUCTION 



Understanding the properties of colloidal systems which are intrinsically out-of- 
equilibrium as they are active, is an important emerging goal in soft matter physics, driven 
both by the potential for modelling motile cells as active colloidal "swimmers", and for 
designing self-propelled particles with novel nanotechnological applications jl-7|. Several 
mechanisms for colloidal self-propulsion have been proposed: these can be broadly grouped 
into those where the colloid exerts force on the surrounding fluid (as for swimming bacteria 
or other motile cells) lj-[5|, and those in which the colloid changes the chemical properties 
of the surrounding medium in an asymmetric way by catalysing a chemical reaction on its 



surface |8 



10| or by secreting some product 
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In this paper, we consider a "chucker": a spherical colloidal particle which produces smaller 
"solute" particles at its surface. Our model colloid is simpler than most cases considered thus 



far 
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-[l5|, since it "chucks" solute particles isotropically, and therefore does not display 
directional motion. However, the chucking drives our system out of equilibrium, giving it 
interesting dynamical properties. We find that the effective diffusion constant of the chucker 
increases with the chucking rate k c for low k c , and eventually decreases for high k c due to 
self-crowding. We further simulate the dragging of a chucker, to obtain a measurement of its 
effective mobility. Comparing this to the effective diffusion constant allows us to show that 
the fluctuation-dissipation theorem does not hold in this non-equilibrium system, and that 
the deviations from the FDT are non-monotonic in the chucking rate. Finally, we consider 
the steady-state positional probability distribution function for a chucker which sediments 
against a planar surface, in the case where the surface is "hard", and in the case where the 
surface absorbs the solute particles. 

Although our model is designed to be generic, our work is inspired by the observa- 



16 



3, 



m some cases 



tion that many bacterial cells secrete extracellular polysaccharides 
in large quantities (for example the bacterium Xanthomonas campestris can produce xan- 
than polymer at a rate as high as 10 4 polymer molecules per cell per second 18|, Il9|). 
These exopolysaccharides can have important effects on the collective properties (e.g. phase 
behaviour) of assemblies of bacteria 20|. Our simulations suggest that in some cases ex- 
opolysaccharide production might also affect the diffusive motion of single bacterial cells. 
We note that the case of a chucker sedimenting onto an absorbing surface might be relevant 
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to exopolysaccharide-producing bacteria located close to a biofilm. 

Several recent works have considered the dynamics of colloids which produce surface 



particles from a theoretical point of view 
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2l| and using simulations jl3|. The key 



concept here is that motion of the colloid can be generated by a concentration gradient of 
solute, produced by the colloid itself. Most previous work has considered the case where 
the colloid produces solute anisotropically across its surface, resulting in directional motion. 
In contrast, we consider the case where solute is produced isotropically, so that the solute 
concentration gradient is on average symmetrical about the colloid. However, as we show 
here, solute production nevertheless has a strong effect on the dynamical fluctuations of 
the colloid. The case of an isotropic chucker has recently been considered theoretically by 
Golestanian fll|. Using scaling theory and linearised hydrodynamics, he predicts inertial and 
diffusive dynamical regimes, separated by an anomalous regime for intermediate timescales 
in which the mean square displacement scales as the 3/2 power of time. He further predicts 
that in the diffusive regime, the effective diffusion constant scales linearly with the chucking 
rate. In this paper, we test these predictions using numerical simulations. 

The question of how a solute concentration gradient influences colloid motion has some 
subtleties. From one point of view, both the colloid and the solute particles are large com- 
pared to the solvent particles; Brownian Dynamics (BD) simulations, in which the solvent 
is replaced by a stochastic force on the particles, are generally believed to work well in this 



regime [22H25| . In such a BD representation, a concentration difference of solute across 
the colloid surface would produce a force on the colloid f = —ksT j iicdS where c is the 
local concentration of solute, the integral is over the colloid surface, n is the normal to the 
surface, ks is the Boltzmann constant and T is the temperature [13j. However, such a BD 
representation does not take account of the details of the solvent motion close to the chucker 
surface. This may be important, since, as emphasised recently by Jiilicher and Prost 21], 



hydrodynamic theory using linearised non-equilibrium thermodynamics predicts that rela- 
tive motion between the colloid and solvent is dependent on the existence of a surface slip 
velocity. In a simple view, the density of the solvent is expected rapidly to compensate any 
imbalances in the solute concentration field to maintain a constant overall isotropic pressure, 
except in a very narrow region around the colloid where molecular interactions between the 
colloid surface and the solute can generate a surface slip velocity. Equivalently, one might 
state that in order to generate a force, the configuration of solute particles relative to the 



3 



colloid surface must change when the colloid moves. Because the solute moves with the 
fluid, such a configurational change can only be achieved if there is a surface slip velocity. 

As discussed by Golestanian et al fill , \\\ . the colloid velocity can be expressed as the 
integral over the colloid surface of a slip velocity v s = — /iVc where c is the solute concen- 
tration, fi the surface mobility fi = kBT\ 2 /r], t] being the viscosity and A a length-scale. In 
the non-permeable case, the latter is the Derjaguin length, describing the range of interac- 
tions between the colloid and solute. The length-scale A is then small, so that the linearised 
theory predicts that chucking will have little effect on colloidal motion. However, Adjari and 
Bocquet J26 1 have shown recently that for solutes smaller than the intrinsic slip length, a 
small degree of surface slip can lead to large enhancements of A. Moreover, for colloids that 
are partially or fully permeable to solvent (osmiophoresis), anything between this result and 



a much larger colloid velocity, equal 
by a mobility factor, is possible LL2 



;o the "naive" kinetic force f = — fcgT J iicdS multiplied 
27J. The motion of a semi-permeable rigid vesicle in 
response to a solute gradient has been analysed in a linearised framework by Anderson ^2S [ . 
In a later work, the applicability of the linearised theory to this type of problem has been 
called into question by an experimental study by Nardi et al ^£7\- These authors measured 
the motion of lipid vesicles in a solute concentration gradient, and found a drift velocity more 
than 3 orders of magnitude faster than that predicted by the linearised theory, suggesting 
that a strong nonlinear coupling exists between osmosis and hydrodynamic flow. 

Our simulations do not include explicit solvent particles. Instead, we model the chucker 
as a diffusing hard sphere and the solute particles as diffusing hard spheres which interact 
only with the chucker. We use a dynamical Monte Carlo (MC) simulation scheme which has 



been shown to be equivalent to Brownian Dynamics for small trial displacements |22M25l. |29|. 
By simulating a passive colloid in a fixed gradient of solute particles, we show that the 
naive kinetic result for the colloid velocity, expected for Brownian Dynamics simulations, 
is recovered in our MC simulations, at least for small solute particles. To estimate the 
likely effects of neglecting the detailed solvent-chucker interactions, we complement our 
MC simulations with a Langevin dynamics approach, in which the chucker is modelled by 
a diffusing point particle which experiences a drift force proportional to the local solute 
concentration gradient. The prefactor linking the force to the density gradient (which is 
proportional to A) can be varied, mimicking the effects that would be expected with explicit 
solvent. The Langevin simulations are in qualitative agreement with our MC results; this 
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agreement is also quantitative, for suitably chosen A, for low chucking rates and when the 
solute particles are small compared to the colloid. 

This work is structured as follows. In the next Section, we describe in detail our simulation 
model and methods. In Section III we report results for a diffusing chucker, a dragged 
chucker and a chucker close to a surface. Finally, Section IV contains a discussion and our 
conclusions. 



II. MODEL AND METHODS 



We model the chucker as a hard sphere of radius Ri, with long-time diffusion coefficient 
D\. With a rate k c , the sphere produces ("chucks") smaller particles (solute) isotropically at 
its surface. The solute particles have radius R 2 and long-time diffusion coefficient D 2 , which 
is related to D\ by Stokes' law: D2/D1 = R\/R 2 . 

In our simulations, the length unit is the chucker radius R\ and we use two different 
values for the solute radius R 2 : R 2 /Ri = 0.1 (R 1 /R 2 = 10) or R 2 /Ri = 0.01 (Ri/R 2 = 100). 



A. Monte Carlo simulations 



Our main simulation results are obtained using a dynamical Monte Carlo algorithm. 
In this algorithm, a particle is chosen at random and subjected to a trial displacement 
S, in a randomly chosen direction and with magnitude uniformly distributed in the range 



— > 0.001-Ri (for the chucker) and — > 0.001i?i \J R\/ R 2 for the solute particles (this 
ratio of displacements is to ensure the correct ratio between the diffusion coefficients D\ 
and D 2 , as given by Stokes' Law). If this trial displacement results in an overlap between 
chucker and solute, it is rejected; otherwise it is accepted. One MC cycle consists of TV 
trial displacements, where iV is the current number of particles in the system. It has been 



shown in previous literature 



22-25 



28 



j|-|30] that this algorithm is equivalent to Brownian 
Dynamics for small trial displacements. An MC cycle can be interpreted as the time unit 
in a dynamical trajectory, where the diffusion constant in reduced length units is given by 
D\jR\ = [MC cycles] -1 . At the beginning of every MC cycle, a new solute particle 

is inserted with uniform probability of k c , randomly at the surface of the chucker. Note 
that because the solutes do not interact, all particle insertions are accepted. All values for 
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the chucking rate k c in our MC simulations are in units of [MC cycles]" 1 . This dynamical 
MC scheme allows us to represent the solute particles as penetrable hard spheres and easily 
to insert new solute particles, both of which would be difficult in a standard Brownian 
Dynamics scheme. 

The chucker and the solute particles interact according to the Asakura-Oosawa-Vrij (or 

33| : the solute particles cannot overlap with the chucker 



penetrable hard sphere) model 
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but do not interact with each other: 








(1) 



if r sc > R1 + R2 



v sc = 00 if r sc <R 1 + R 2 



where v ss and v sc are the solute-solute and solute-chucker interaction energies, and r sc is the 
distance between the centres of the solute and chucker particles. 

Figure [T] shows snapshots from our MC simulations. To prevent the number of solute 
particles from growing indefinitely during the simulations, solutes need to be removed at 
some large distance from the chucker. Our system reaches a steady state when the rate of 
solute particle generation at the chucker surface equals the rate of solute absorption at the 
boundary. We achieve this in two different ways. For simulations with a fixed wall (as shown 
in the right-hand panel of Figured]), we define a fixed cubic simulation box and remove any 
solute particles which stray outside this box. For simulations in the bulk (as shown in the left 
panel of Figured]), we remove any solute particles that stray outside a cubic box with sides 
of length 2L, centred on the current position of the chucker. The box size is large enough 
to have no effect on the results reported here (L = QR\). We have also verified that for our 
parameter set, the results are independent of the box size. Figure [2] shows the steady-state 
concentration profile of solute particles around the (moving) chucker, for several different 
values of L: for L > 6R±, the solute concentration profile is virtually independent of L. 

B. Langevin Dynamics simulations 

To facilitate the simulation of larger systems, and to better understand the underlying 
physics, we have also simulated our system using an overdamped Langevin dynamics algo- 
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rithm 



34), in which the chucker is represented as a point particle and the solute particles 



are represented as a space and time-dependent density field, the gradient of which provides 
a force on the chucker. In a mean-field approach, we assume that the density of solute parti- 
cles can be obtained from the time-dependent analytical solution of the free-space diffusion 
equation with a point source (the chucker), integrated over the previous trajectory of the 
chucker. This approach assumes that the thermal fluctuations of the chucker are slower 
with respect to the solute diffusion, and further neglects thermal fluctuations in the solute 
density. A similar approach has recently been used to simulate a microorganism performing 
chemotaxis by Sengupta et al [3|- 

In our Langevin dynamics simulations, the force f on the "point chucker" depends on the 
concentration profile c(r) of the solute according to 
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f = -AVc(r) (2) 

where the prefactor A is determined by the physical process by which an osmotic pressure 
gradient is converted into colloidal motion. This force produces a drift velocity v which 
is related to the force f by f = ksTv/Di = 6irr]Ri\. By tuning the prefactor A, we can 
tune the drift velocity v. If the solvent- chucker interactions were to be properly taken into 
account, linearised theory s ugge sts that the drift velocity would be given by Golestanian's 
relation v s = — /csTA 2 Vc/?7 where (as discussed earlier) A is related to the lengthscale 
of solvent-colloid interactions. By combining this relation with Eq.(T2]), we can obtain an 
"effective" value of A for our Langevin dynamics simulations: A = \F~Aj {QuRiksT). A full 
treatment, including solvent, with this value of A, would give the same drift velocity as 
our Langevin dynamics simulations with prefactor A. By varying the prefactor A, we can 
vary the effective A value. By choosing A values corresponding to realistic solvent-colloid 
interaction lengths, we should obtain results closer to those that would be expected for a 
full treatment including solvent. In contrast, in our Monte Carlo simulations, the value of 
A is fixed. Here, the net force on the chucker produced by the imbalance of solute collisions 
across the colloid surface 13j corresponds to a value A = (4/3)7ri?f &bT, or alternatively an 



effective A = 

The concentration profile c(r) is not represented explicitly in our simulations but is cal- 
culated analytically from the time history of the position of the point chucker: 

k c [ dt'G(r-r',t-t'), (3) 
Jo 



c r 



where r = r(t), r' = r(t') and G(r — r', t — t') is the Green's function for the diffusion 
equation in free space with a sink at infinity (valid for comparison to our MC simulations 
when the MC boundary position L is large): 



G(r-r',t-0 = 8(7rjD2( /_ f))3/2 exp 



(x - x') 2 + (y - y') 2 + (z - z'f 



4£> 2 (t-f) 

where D 2 is the diffusion constant for the solute. We now define F(r — r',t — t') = 
VG(r — r',t — t'). Since in our Langevin dynamics algorithm the particle moves in dis- 
crete steps, we can replace the integral in Eq.([3]) by a sum of F(r — r',t — t') over the 
previous trajectory of the particle, giving us a history-dependent expression for the force on 
the particle: 

M 

f it) = -Ak c AtF (r(t) - r(t - iAt), iAt) (5) 
i=i 

where At is the timestep. The sum in Eq.(j3]) is carried out only over the previous M steps, 
where MAt should be longer than the timescale over which the solute diffuses away. The 
equation of motion for the point chucker is then: 

r(t + At) = r(f) + ^Atf (f) + £ (6) 

where D\ is the diffusion constant of the chucker and £ represents Gaussian white noise with 
zero mean and variance 2D\At. Care is needed when selecting the parameters D±, D 2 and k c 
for comparability of the results of the Langevin and MC simulations, because of the different 
time units of the two algorithms; however, dimensionless combinations of parameters can be 
defined which are equivalent between the two simulation schemes. 

In our MC simulations, we expect that at high chucking rates, the solute particles will 
tend to form a cage around the chucker, restricting its diffusion. This effect is not included 
in the Langevin dynamics simulations where we consider a point chucker and a continuous 
solute field. The Langevin Dynamics simulations are therefore only valid for low chucking 
rates, where the solute concentration is small. 



III. RESULTS 

In Sections IIII Al and IIII B\ we first discuss the behaviour of a chucker in free space, with 
no external applied forces. We then determine the response of the chucker to a dragging 
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force in Section IIII C\ and finally we discuss the behaviour of a chucker which sediments 
against an absorbing or hard-planar surface in Section IIII Dl 

A. Solute concentration profiles 

We begin by investigating the concentration profile of solute particles around the chucker 
in our MC simulations when R1/R2 = 10, for a fixed chucker. Since the solute particles do 
not interact with each other, we expect the steady-state solute concentration profile to cor- 
respond to the solution of the diffusion equation ([3]), using the Green's function appropriate 
for a sink at distance L from the chucker. This solution (for the solute concentration c) is 
c ( r ) = 4wb~ (r ~~ r)> so solute volume fraction is given by 0(r) = (}■ ~ t)- 

We note that this is not a true volume fraction, since solute particles do not interact and can 
overlap in our model (so that <fi can be greater than one); but it does give some indication 
of the density of the solutes, for comparison to real systems. 

Figure [3] shows a comparison between the radial profiles for <fi obtained from our simu- 
lations with the analytical expression, for L = AR\. The simulation data agrees very well 
with the analytical expression; the small discrepancy close to the boundary (r/Ri = 4) 
occurs because the analytical expression is computed for a spherical boundary whereas the 
simulation has a cubic boundary. 

We are also interested in how closely the solute particles are packed around the chucker: 
if the volume fraction close to the chucker surface is too high then the the Asakura-Oosawa 
approximation fails to accurately describe a real solute (e.g. polymer coils), since it neglects 
solute-solute interactions. Figure H] shows the volume fraction of solute particles in spher- 
ical shells centred on the chucker in our MC simulations, computed for a freely diffusing 
chucker with R1/R2 = 10, over the full range of chucking rates used in our simulations 
(0 < k c R\jD\ < 600). The maximum volume fraction obtained close to the chucker is 
about 0.12, and for most chucking rates it is much lower. This suggests that our neglect of 
solute-solute interactions is justified over the range of parameters used in this work. 
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B. Chucker diffusion 



We now discuss the diffusive behaviour of an isolated chucker, using MC simulations 
with two different ratios of the chucker and solute radii: R1/R2 = 10 and R1/R2 = 100, for 
various chucking rates k c . In our simulations, the solute particles are produced isotropically 
so there is no net directional motion. Rather, we observe diffusive behaviour, in which the 
mean square displacement (r(t) 2 ) is linear in time. This allows us to compute an effective 
diffusion constant D e g by linear fitting of the long-time behaviour of (r(t) 2 ). 

Figure [5] shows that in our simulations the effective diffusion constant _D e ff of the chucker 
varies non-monotonically with the chucking rate k c . Diffusive motion of the chucker is 
enhanced as the chucking rate increases, for low chucking rates. At high chucking rates, 
diffusive motion is inhibited due to crowding by the solute particles, so that D e g decreases 
with k c . We plot the dimensionless quantity D e s/Di, where Di is the diffusion constant of 
the chucker for zero chucking rate, versus k c R\jD\. The latter turns out to be the most 
relevant dimensionless parameter for comparison to theory: see Ref [ll| and Eq.((7]). R\jD\ 
is the typical time for the chucker to diffuse the radius of a solute particle; this ratio therefore 
measures the relative magnitudes of this time and the inverse of the chucking rate. 

We believe that the mechanism behind the increase in D c g with chucking rate, for small k c , 
is as follows. A stationary chucker would be surrounded by a symmetric, concentration field 
of solute which decreases with distance from the chucker. A small random displacement of 
the chucker, however, positions it asymmetrically in the solute concentration field. Osmotic 
effects will then tend to cause the chucker to move further in the same direction, leading to 
a non-thermal osmotic enhancement of the effective diffusion constant of the chucker. 

Comparing the results in Figure for the two different ratios of chucker and solute 
radius, we find that the small k c regime, where _D e ff/-Di increases with k c R\jT)\, appears 
not to depend strongly on the size of the solute - although the inset, plotted on a linear 
scale, suggests that the two plots do not actually collapse convincingly onto one curve. On 
the other hand, our data suggests that the regime, for large k c , in which diffusion decays 
with chucking rate, does depend strongly on the size of the crowder: crowding effects kick 
in at smaller k c R\jD\ for the smaller solute particles. This suggests that the dependence of 
the crowding effect on solute radius is less steep than ~ R\. 

Golestanian's recent theoretical work which uses linearised hydrodynamics, predicts 
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that Defi/Di should be given by: 



D 



eft 



A\ 4 /3 



1.17810 I — , 

R\) \7T 



(7) 



where we have substituted the Stokes-Einstein relation {D_\ = 6 fc ^ T ) and the relation 



6-kRiT] > 

/i = ksTX 2 /?] into the expression given in Eq.(7) of Ref ll|. In our MC simulations, 
the imbalance of solute collisions across the chucker surface is expected to produce a net 
force which corresponds to a prefactor A = (A/3)iiR\kBT in Eq. (12|) . This will produce a 
drift velocity v = f / (QirqRi) . Linearized theory, which accounts for solvent- colloid surface 



interactions, predicts a colloid velocity v s = — A;^TA 2 Vc/?7 11] . Matching these two rela- 
tions, we find that our MC simulations have an effective A value A = \/2Ri/3. Substituting 
this into Eq.(17j), we obtain a prediction for D e Q in our MC simulations. 

The inset to Figure [5] shows a comparison between this prediction and our MC data, for 
small chucking rates. For R1/R2 = 100, we obtain an excellent fit, whereas for the larger 
solute particles, R1/R2 = 10, the fit is much less convincing. We will revisit the effect of 
solute size in Figure [7J We note that Golestanian jll| also predicts a regime in which the 
mean square displacement (r 2 ) is proportional to t 3 / 2 , for timescales between the inertial and 
diffusive regimes, for an isotropic chucker. We do not observe this regime in our simulations: 
for ratios of R1/R2 > 5 the mean square displacement is linear in time. 

To further investigate the relationship between chucker diffusion and solute concentration 
gradient, we performed additional MC simulations for a passive colloid, of radius R\, in a 
constant density gradient of Asakura-Oosawa solute particles (FigEJ). In order to maintain a 
constant density in steady state, we introduced a planar source of solute particles (of radius 
R2) at the top of the simulation box, and a planar absorbing sink at the bottom. Solute 
particles are randomly seeded on the source plane at a constant rate k c . 



The analysis in Ref. 12j can be generalised to yield the following prediction for the 
steady state velocity of a colloidal particle in a constant concentration gradient (neglecting 
hydrodynamic effects): 

A-KR 2 R\k c 

where the top and bottom planes have a surface area of L 2 , and we have again used the 
expectation that A = (y/2/3)Ri in our simulations. 

Figure [7J shows a comparison between the analytical prediction (jSJ) and our simulation 
results, for R1/R2 = 10 (panel a) and R1/R2 = 100 (panel b). As in the inset to Figure 0, 
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the agreement with the theoretical prediction is much better for the smaller solute radius. 
This observation is not unexpected, given that the smaller solute radius is closer to the 
continuum limit. What is perhaps surprising is the extent of the deviation from the theory 
when R1/R2 = 10. For this solute radius, our simulation results deviate from the analytical 
prediction by about an order of magnitude. 

An important drawback of our MC simulation scheme is that it does not account for the 
details of the solvent interactions with the chucker surface. As discussed in Section [TJ these 
effects may drastically reduce the magnitude of solute-induced effects on colloid mobility. In 
our MC simulations, an imbalance of solute concentration across the chucker surface leads to 
fewer overlaps ("collisions") on one side of the chucker than the other, and hence net motion 
down the concentration gradient. We expect this to correspond to Eq.([2]) with a prefactor of 
A = (4/3)irRlkBT (simply summing collisions across the colloid surface). We would expect 
hydrodynamic effects at the chucker surface (if included) to decrease this prefactor: for 
a solvent-impermeable chucker with no slip boundary conditions, linearised hydrodynamic 



theory predicts that the prefactor decreases by a 



for chucker- solute interactions! is ~ Ro 
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actor ~ \ 2 / R\ where A [the lengthscale 



30|. If, on the other hand, the chucker is 



partially or totally permeable to solvent, the decrease in the prefactor due to these effects 
should be less severe 
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27|. 



As discussed in Section III B[ we can model the effects of solvent-chucker surface interac- 
tions approximately using Langevin dynamics simulations in which we vary the magnitude 
of the prefactor A: this corresponds to varying the hydrodynamic lengthscale A. Figure [8] 
shows the effective diffusion constant D c s, computed as a function of k c using Langevin 
dynamics simulations with different prefactors, for Ri/R 2 = 10. The same qualitative be- 
haviour is observed as in our MC simulations (compare to Figure [5]): again, D c g/Di increases 
with k c for small chucking rates. For large chucking rates, the effective diffusion constant de- 
creases in the MC simulations due to crowding effects which are not included in the Langevin 
Dynamics simulations 

For small k c R\j D\, we can make a quantitative comparison between the MC and Langevin 
simulation results, as shown in Figure [9j Here, the Langevin simulations were run with the 
"naive kinetic" prefactor A = (4/3)-7ri?f/c^T. Excellent agreement is obtained between the 
two types of simulation. This reinforces our earlier conclusions that the MC simulations 
for small R2 behave as would be predicted by hydrodynamic theory with a hydrodynamic 
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lengthscale A = (a/2/3) Ri- 



C. A dragged chucker 

For passive systems, such as colloidal particles, the fluctuation-dissipation theorem (FDT) 
states that the mobility coefficient V is fundamentally related to the diffusion constant D by 
T = D/ksT. For active systems, there is no reason why the FDT should hold. Deviations 
from the FDT not only serve to demonstrate that the system is out of equilibrium, but can 
also yield insight into the underlying physics. 

To investigate the relationship between mobility and diffusion for the colloidal chucker, 
we carried out MC simulations in which a chucker (with R1/R2 = 10) is dragged in a fixed 
direction by a constant force F (for example, using optical tweezers). We measure the 
velocity v of the chucker in the direction of the force. The mobility coefficient T is defined 
by v = TF; this can be compared with the value of -D e ff obtained in Figure [5j In this work, 
we use values of F which approximately correspond to the magnitude of the gravitational 
force on a bacterium. 

Figure [TUl shows the steady state chucker velocity v, as a function of the chucking rate, 
expressed as the dimensionless ratio k c R\jD\, for two different values of the pulling force F . 
This velocity depends strongly and nonmonotonically on the chucking rate, in a markedly 
different way to the effective diffusion constant of Figure [5j 

Combining the results of Figures [5] and [101 we obtain Figure HH which describes how the 
ratio ThsT / D c r depends on k c R\/D\. Deviation of TksT / D c r from unity (green dashed 
line) corresponds to violation of the FDT. We observe firstly that deviations from FDT in 
this system are strong (up to a factor of 3) and nonmonotonic. Depending on the chucking 
rate and the pulling force, mobility may dominate diffusion, or vice versa. As might be 
expected, the deviations from FDT are stronger for the larger pulling force. While it is 
not surprising that the FDT does not hold in this nonequilibrium system, what is perhaps 
unexpected is that the relationship between the mobility and diffusion constants has such a 
complex dependence on both the chucking rate and the applied force for this rather simple 
model. 

The trends observed in Figure [TT] can be tentatively related to changes in the solute 
configuration around the chucker, for different chucking rates and pulling forces. Naively, one 
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might imagine that for small chucking rates, where the density of solutes is low, the pulling 
force might tend to displace the chucker relative to its symmetric concentration gradient 
of solute, so that more solute would be behind it than in front, resulting in an osmotic 
enhancement of the mobility. This is what one would expect from an advection-diffusion 
equation with a source moving at a finite Peclet number (see e.g. Ref. ^6|). However, for 
large chucking rates, the density of solutes near the chucker is higher, and one might imagine 
that the pulling force would lead to an increased collisions of the chucker with nearby solute 
particles, decreasing the mobility. 

To test this hypothesis, we plot in Figure [12] the "volume fraction", as defined by = 
(4/3)7ri?2C, of solute particles within the conical volume supported by the angle 9 = tt/9, 
in front of and behind the chucker (the cone is extended all the way to the edge of the 
box). We first note that there is a difference in the volume fraction of solute in front of and 
behind the chucker, for large enough chucking rate. These effects are more pronounced for 
the larger pulling force. For F = 7 .OksT / Ri, the solute concentration is larger behind than 
in front of the chucker, for chucking rates in the range 50 < k c R\jD\ < 200. We would 
expect this to cause an enhancement of mobility over diffusion; however, this range of k c 
values is beyond the peak in r/cgT '/ D e g (as shown by the green dashed line). For higher 
chucking rates, we would expect to see a larger volume fraction in front of the chucker, if the 
solute particles are inhibiting chucker motility. There is some indication that this happens 
for the highest chucking rates (not shown in Figure [T2l . It is possible that this crowding 
effect plays a more important role for the smaller value of the force F = 3.0&bT/.Ri, where 
the pulling force is not strong enough to overcome the crowding; for the stronger pulling 
force F = 7 tiksT / R\, we speculate that the pulling force is strong enough to overcome this 
effect. Further work will be needed to fully understand the data in Figure [TTJ it is clear, 
however, that the nontrivial relationship between mobility and diffusion in this system is 
caused by a complex interplay between crowding and osmotic effects. 



D. A sedimenting chucker 

We now consider the case of a chucker which experiences a constant force F perpendic- 
ular to a planar surface which either absorbs the solute particles or behaves as a hard wall. 
This models the gravity- driven sedimentation of a chucker onto a surface, in the case where 
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gravitational effects on the solute particles are negligible. Although this model is highly 
simplified, the sedimentation of exopo 



ysaccharide-producing bacteria onto surfaces is im- 
portant in the formation of biofilms [37|. Here, the hard surface case might mimic the initial 
stages of biofilm formation, when bacteria sediment onto an uncolonised surface, while the 
absorbing surface might correspond to a situation where bacteria sediment onto an already 
formed biofilm which tends to absorb exopolysaccharide. 

We determine the steady-state probability distribution function for the chucker position 
relative to the surface in our MC simulations. This would correspond to an experimental 
sedimentation density profile, in the case of low chucker density where interactions between 
chuckers can be ignored. In all our simulations, we take R1/R2 = 10, F = 7.0/csT/i?i and 
initialise the simulations with the centre of the colloid at a distance z = I.6-R1 from the 
surface; we then begin to observe the system once a steady state has been reached (i.e. after 
9 x 10 5 Monte Carlo cycles). This choice of F approximately corresponds to the gravitational 
force experienced by a typical bacterium. 

Figure [13] shows the steady-state probability distribution P(z) for the position of the 
centre of the chucker relative to the surface, which is at z = 0. In Figure ITBT a). the wall acts 
as a sink for solute particles. This results in a buildup of the chucker probability density 
close to the surface: the higher the chucking rate, the more tightly the chucker tends to 
approach the surface. This effect can be understood in terms of an effective attractive 
interaction between the chucker and the surface, in which local depletion of solute near 
the surface causes an osmotic imbalance which squeezes the chucker against the surface. 
Figure fl37 b) shows the corresponding probability density profiles when the planar surface 
instead acts as a hard wall. Here we find two regimes, depending on the chucking rate. For 
small chucking rates (k c R\/Di < 1), the probability density decreases monotonically away 
from the surface, with a steepness that decreases as k c increases (for k c = we obtain a 
Boltzmann distribution). However, as k c increases further, the profile develops a minimum 
close to the surface and a peak at some distance from the surface. This peak recedes from 
the surface with increasing k c . For the hard surface, there is a solute-mediated osmotic 
repulsion between the chucker and the surface. The crossover between the two types of 
probability density profile occurs when this repulsion becomes strong enough to overcome 
the gravitational force pulling the chucker towards the surface. In the high k c regime, the 
peak in P(z) occurs where the gravitational and osmotic repulsion forces balance. We note 
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that for the hard wall, the solute density profile, and hence the steady-state chucker location, 
depends on the position of the absorbing box boundary, since a linear concentration gradient 
of solute arises between the colloid and the box boundary. However, the qualitative features 
of our results are independent of the boundary position. 

IV. DISCUSSION 

In this work, we have used computer simulations to study the physics of an active colloid 
(a "chucker") of radius Ri, which produces solute particles of radius R2 isotropically on 
its surface, at a rate k c . Most of our simulations have used a Monte Carlo scheme in 
which the chucker is modelled as a hard sphere and the solute particles as Asakura-Oosawa 
(i.e. mutually non-interacting) spheres. Our MC simulations were complemented by an 
overdamped Langevin dynamics approach in which the chucker is represented as a point 
particle which experiences a force proportional to the gradient of a concentration field of 
solute, this field being computed analytically using the Green's function for free diffusion. 
The Langevin dynamics approach has the advantage that we can tune the prefactor that 
links the osmotic force to the solute concentration gradient: this provides a crude way to 
model hydrodynamic effects close to the chucker surface. 

Our MC simulations of an isolated chucker in free space show diffusive behaviour, from 
which we extract an effective diffusion constant D c g. This varies nonmonotonically with the 
chucking rate. For small values of k c , solute chucking enhances diffusion. We interpret this 
as being due to the chucker being "pushed along" by a self-generated solute concentration 
gradient; the same behaviour is recovered in our Langevin dynamics simulations. An increase 
in -DcfT as a function of k c was recently predicted analytically by Golestanian . In our MC 
simulations, when the chucking rate becomes large, the chucker tends to become caged by 
newly generated solute particles, hampering its diffusion. The result is therefore a peak in 
D e g at intermediate chucking rate. In our simulations, the solute particles do not interact. 
For more realistic, interacting solutes, we might expect crowding to kick in at lower chucking 
rates, shifting the peak in D e g towards smaller k c . 

A quantitative comparison between our MC simulation results and Golestanian's theoret- 
ical prediction reveals that the simulations are in agreement with theory for a hydrodynamic 
length A = (v2/3)i?i, for large values of the ratio between the sizes of the chucker and of the 
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solute particles, R1/R2 = 100; agreement is less convincing for R\jR 2 = 10. This is the value 
of A we would predict by matching the drift velocity predicted in a "naive" osmotic view, in 
which the colloid moves because of an imbalance in the number of solute collisions across its 
surface, with the colloid velocity, expressed as a function of A, predicted by hydrodynamic 
theory. When we simulate a passive colloid in an externally imposed concentration gradient 
of solute particles, we also find good agreement with the hydrodynamic theory prediction 
with A = (V2/3)i2i for Ri/R 2 = 100, but for R\jR 2 = 10, the diffusiophoretic mobility is 
smaller than might naively be predicted from osmotic pressure arguments. We speculate 
that this poor agreement for large solutes may be due to local fluctuations of the solute con- 
centration field, neglected in the theory. We can also make a quantitative comparison, for 
small chucking rates, between our MC and Langevin simulation results. For R\/R 2 = 100, 
the two simulation schemes produce effective diffusion constants in good agreement, for a 
choice of Langevin prefactor of A = (4/3)71"/?^ /cgT, corresponding to A = (v / 2/3)i?i. 

One of the initial motivation for this study was to investigate possible effects of ex- 
opolysaccharide production on bacterial diffusivity. The bacterium X. campestris produces 
xanthan polymer at a rate as hi gh as 10 4 polymer coils per cell per second (as estimated 



from the data given in Refs. 
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19[). Assuming a ratio R\/R 2 = 10 — 100, corresponding 



to a polymer gyration radius of 10-100 nm, hence to a value of k c R\jD\ ~ 1 — 100, and 
referring to Fig. we find that the effective diffusion constant of the chucker should be 
significantly different from the infinite dilution value at this chucking rate. On the other 
hand, the bacterium Sinorhizobium meliloti, a less prolific^ producer of exopolysaccharide, 

, leading to a value of k c R\jD\ 



produces less than 1 polymer per bacterium per second 
of 10 -4 — 10~ 2 (assuming the same dimensions for the polymer), so that for this bacterium, 
we expect the effect of polymer production on diffusion to be essentially perturbative. It 
would be interesting to compare the diffusion constants of mutants of X. campestris which 
are altered in their exopolysaccharide production rates: our simulations suggest a change in 
-D e ff of up to about 50% if exopolysaccharide production were to be eliminated. One would 
need to be careful in such an experiment to avoid the buildup of polymer in the bacterial 
suspension, leading to crowding and a subsequent decrease in diffusion. 

The chucker is an example of a nonequilibrium, "active" system, in which statistical 
mechanical relations such as the Fluctuation-Dissipation theorem (FDT) may not hold. 
Comparing the effective diffusion constant D\ to a mobility coefficient computed by dragging 
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a chucker with a constant force in our MC simulations, we find that indeed this system 
violates the FDT strongly, in a highly nontrivial way. The complex relationship between 
mobility and diffusion in this system is due to an interplay between crowding and osmotic 
effects which depends sensitively on the applied force and chucking rate. 

Finally, we have studied the behaviour of a chucker which sediments onto a planar surface 
which either absorbs the solute particles or behaves as a hard wall. When the surface is 
absorbing, a solute-mediated effective attraction (which increases with increasing k c ) pushes 
the chucker towards the surface. However, when we consider a hard surface, a solute- 
mediated effective repulsion exists. For small k c this simply results in a decrease in the 
gradient of the chucker probability density profile near the surface. For large k c , the repulsive 
interaction dominates the gravitational force near the surface, resulting in a peak in the 
probability density for the chucker position at some distance from the surface. 

This work was motivated by recent experimental, theoretical and simulation work on 
colloidal "swimmers" which catalyse a chemical reaction, or secrete some product, anisotrop- 
ically across their surface js i5|, by recent theoretical work on isotropic chuckers and 



by the biological example of exopolysaccharide-secreting bacteria 
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17j. However, the 



physical mechanisms explored here also have wider relevance in soft matter physics. As an 



example, Jiang et al 



39| recently reported a system where migration of colloidal particles 



was mediated by inhomogeneities in polymer concentration, which were in turn induced by a 
temperature gradient. Dynamic depletion interactions between colloidal particles, in which 
one of the components of a phase separating binary fluid wets the colloid surface (creating 
a local depletion near the colloid), predicted in simulations by Araki and Tanaka ^40 1 . are 
also closely related to the phenomena studied here. 

We hope that the results reported here will make a useful contribution to ongoing dis- 
cussions concerning the effects of osmotic gradients on colloidal dynamics, as well as link- 
ing these discussions to the physics of biological systems (particularly exopolysaccharide- 
producing bacteria) and other nonequilibrium soft matter systems. 
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Figure 1: Snapshot of a freely diffusing chucker (left), and a chucker close to a hard-planar surface 
(right). 
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Figure 2: Volume fraction (<ft = (A/fyitR^c (where c is concentration) )of solute particles surrounding 
a motile-chucker (R1/R2 = 10) with k c R\/D2 = 600 (k c = O.OOlfA/Ccyc/es]" 1 ) and various box 
size L. 
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Figure 3: Comparison between analytical expression (green dashed-line) and simulated data (solid 
black line) for the "volume fraction" <fi = (4/3)7ri?2 c (where c is concentration) of solute particles 
surrounding a non-motile chucker (R1/R2 = 10) with different values of k c (from top to bottom, 
k c = 0.00001, 0.0001, 0.001 [MC cycles] -1 ). The corresponding values of the dimensionless ratio 
k c Rf/D2 are given in the legend, for comparison with the analytical expression. Note that since 
the solute particles can overlap, <j> is not a true packing fraction. 
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Figure 4: "Volume fraction" <\> = (4/3)-7ri?2C (where c is concentration) of solute particles surrounding 
the chucker (R1/R2 = 10) as a function of the distance from the centre of the chucker (expressed in 
units of Ri). The corresponding values of the dimensionless ratio k c R\/D2 are given in the legend 
(note that the solute particles can overlap). 
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Figure 5: Effective diffusion coefficient D e g of a chucker (normalised by the diffusion constant in the 
absence of chucking D\) as a function of the dimensionless ratio k c R\/D\, for R1/R2 = 10 (black 
circles) and R1/R2 = 100 (red squares). Note the logarithmic scale on the horizontal axis. The 
inset (on a linear scale) shows a comparison between the same data for small k c and the theoretical 
prediction, Eq.([7]), with A = \/2i?i/3 (green dashed line). 
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Figure 6: Snapshot of a passive colloid in a concentration gradient of Asakura-Oosawa solute 
particles: solute particles are generated at the top planar source, and removed at the bottom 
planar sink (the positions of these planes are indicated by the blue circles). 
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Figure 7: Sedimentation velocity for a passive colloid in a concentration gradient of solute, generated 
by a planar source and absorbed at a planar sink. The dimensionless combination k c R\jD\ measures 
the solute concentration gradient (note the colloid does not chuck in these simulations; solute is 
introduced into the simulation box at the source plane with rate k c ). The black circles indicate 
analytical results from Eq.([5]), whereas the red squares are simulation data, (a): the case R1/R2 = 
10 (b): the case R 1 /R 2 = 100. 
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Figure 8: D c q computed with Langevin dynamics simulations. The Langevin dynamics simulations 
are performed with various prefactors A = aksTRf, where /c#T = 1, and the values of a are given 
in the legend. 
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Figure 9: D c r/D\ versus the dimensionless ratio k c R\jD\ for a chucker with R1/R2 = 100: the 
black continous line represents the MC data and the red dashed line the Langevin results with 
A = (4/3)71"^ /cbT. The plot shows data only for small values of k c R\jD\. 
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Figure 10: Asymptotic velocity v versus the dimensionless ratio k c R\jD\ for a chucker under an 
external force of F = 3.0 k B T / R x or F = 7.0 k B T / R x (Ri/R 2 = 10). 
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Figure 11: T ^ bT versus the dimensionless ratio k c R\lD\, for a chucker with R1/R2 = 10. D e g is 
measured from the mean square displacement of a freely diffusing chucker (Figure [5]) whereas T is 
obtained by dragging the chucker with an external force of F = 3.0ksT/Ri or F = 7.0ksT/Ri 
(Figure [10]) . The dashed line indicates the FDT prediction for a passive system. 
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Figure 12: "Volume fraction" <p = (4/3)7ri?2 c (where c is concentration) of solute within a cone of 
angle = tt/9 at the front (red crosses) and at the back (black dots) of the pulled chucker, as a 
function of the dimensionless ratio kcR^/Di, for a chucker with R1/R2 = 10 pulled by an external 
force of F = S.OksT/Ri or F = 7 '.OksT / R\. The green dotted line represents the position of the 
peak of TksT / D e R in Figure [TTJ For F = 3.0ksT/ R±, the peak of TkBT/D e g is around 360 and 
not shown. Note also that the solute particles can overlap. 



32 



k.R, 2 /D 


=0.0 


k.R, 2 /D 


=0.06 


- k R 2 2 /D 


=0.6 


k,R, 2 /D 


=6 




Figure 13: Steady-state probability distribution for the z-coordinate of the centre of the sedimentin^ 
chucker, for R1/R2 = 10 and various chucking rates (indicated in the legend). 
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